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Abstract. A new algorithm for real root isolation of polynomial equa- 
tions based on hybrid computation is presented in this paper. Firstly, 
the approximate (complex) zeros of the given polynomial equations are 
obtained via homotopy continuation method. Then, for each approxi- 
mate zero, an initial box relying on Kantorovich theorem is constructed, 
which contains the corresponding accurate zero. Finally, Krawczyk in- 
terval iteration with interval arithmetic is applied to the initial boxes 
so as to check whether or not the corresponding approximate zeros are 
real and to obtain the real root isolation boxes. Meanwhile, an empirical 
construction of initial box is provided for higher performance. And our 
experiments on many benchmarks show that the new hybrid method is 
more efficient, compared with the traditional symbolic approaches. 



1 Introduction 

The Real Roots Isolation of polynomial equations is a procedure that uses disjoint 
regions to isolate all the distinct real roots of polynomial equations, with only one 
root in each region. Formally speaking, let F — (/i, /2, . . . , f n ) T be polynomial 
equations defined on _R™, i.e. fa € R[xi, X2, • ■ • , x n ]. Suppose F(x) = has only 
finite many real roots, say C , C^ 2 : • • • > £ ■ The target of real root isolation is 
to compute a family of regions Sx,S2,-.-, S m , Sj C R n (l < j < m), such that 
fG) G Sj and Si n Sj = 0(1 < i,j < to). Usually, we use rectangular boxes to 
denote the regions above. So we often call these isolated boxes intervals in this 
paper. 

Real root isolation is an important problem in symbolic computation. It can 
be viewed as a kind of exact algorithm for solving equations since no root formula 
is available in general situation. It is also a critical part of some other important 
algorithms, such as CAD and real root classification for semi-algebraic system, 
etc. Improvement on real root isolation will benefit all of these algorithms. 

We impose some hypothesis on the problem discussed here. First is that the 
system is square, i.e. the number of equations is the same as that of variables. 
Then we only handle the systems with finite many roots. Positive dimensional 
solution is beyond the scope of this paper. Moreover, we suppose that the Jaco- 
bian matrix of F is non-singular at each root of F(x) = 0. So we only deal with 



the simple root cases. For the singular situation, the deflation method [9,10,21,7] 
can be applied, which is one of our ongoing work. 

Most of the previous real root isolation algorithms are based on symbolic 
computations. For instance, the Uspensky algorithm [6] based on Descartes' 
rule is for polynomials in one variable. In multi-variable scenario, we have "First 
algorithm" based on monotonicity [24] and "Second algorithm" based on "upper- 
lower bound" polynomial [23]. There are also some other algorithms based on 
different techniques, see for example [5,3,13,4]. 

An advantage of those symbolic methods is that exact results can be obtained 
since they use symbolic computation and some of them can be extended to semi- 
algebraic systems. However, there are also some disadvantages. Some of these 
method could only handle the isolation of complex zeros. And some of them need 
to triangularize the system first, which is unacceptable in computation when 
the system is complicated sometimes, such as more variables or high degrees. 
While some methods that do not use triangularization have to give a huge initial 
interval to include all the real roots [25,26], which is extremely inefficient. 

In order to avoid these problems and design a new algorithm that could 
efficiently solve more complicated systems and provide accurate interval results, 
we employ hybrid computation to take both the advantages of symbolic and 
numerical methods. 

The basic idea of this paper is to use numerical method to obtain all the ap- 
proximate zeros of polynomial systems, including possible non-real ones. With 
these approximations, small initial intervals which contains the corresponding 
real roots are constructed. We then apply symbolic method to these initial in- 
tervals to verify whether there is a real root in it or not. The main method we use 
in numerical computation is homotopy continuation, and for symbolic process 
we use Krawczyk iteration. 

Most of the work in this paper comes from [16]. In section 2, we will introduce 
some preliminaries, including homotopy continuation and interval arithmetic. A 
new real root isolation algorithm is discussed in section 3. To test our new 
method, our experimental results on benchmarks together with comparison and 
analysis will be presented in section 4. Finally, there is a summary in section 5 
and some future work will also be discussed. 

2 Preliminary 

We introduce in this section some basic theories and tools that would be used 
in our algorithm. 

2.1 Homotopy Continuation Method 

Homotopy continuation method is an important numerical computation method, 
which is used in various fields. We only treat it as an "algorithm black box" here, 
where the input is a polynomial system, and the output is its approximate zeros. 
Please find the details about the theory in [11,17]. 



For our purpose, it is convenient to apply some existing softwares, such as 
Hom4ps-2.0 [12], PHCpack [18] and HomLab [19]. 

In our implementation, we use Hom4ps-2.0, which could return all the ap- 
proximate complex zeros of a given polynomial system efficiently, along with 
residues and condition numbers. 

2.2 Interval arithmetic 

Interval arithmetic plays an important role in real root isolation algorithms 
[25,26,23]. The two main differences between our new algorithm and the tradi- 
tional ones in [25,26] arc: 1) Verification only carry out on the localized "small" 
intervals; 2) symbolic computation is replaced with floating point numerical com- 
putation. 

Most of the interval operations in this paper's algorithms are based on 
Rump's floating point verification work [15] and accomplished by using the Mat- 
lab package Intlab [14], including interval arithmetic operations and Jacobian 
matrix, Hessian matrix calculations. 1 So we will focus on the real root isolation 
algorithm. 

Basic concepts We introduce some basic interval arithmetic theories in this 
section. Sec reference [20] for more details. 

For given numbers x, x £ R, if x < x, then we call 

X = [x,x] = {x £ R\x < x < x} 

a bounded close interval, or interval for short. Denote by I(R) the set of all the 
bounded close intervals on R, and 1(A) = {X £ I(R)\X C A} all the intervals 
on A C R. Especially, if x = x, we call X a point interval. 
For intervals, there are some common quantities: 

midpoint mid(X) = (x + x)/2 
width W(X) =x-x 
radius rad(X) = ±W(X) 
low end point inf (X) = x 
high end point sup(X) = x 

Obviously we have X = [mid(X) — rad(X), mid(X) + r&d(X)]. An interval is 
usually expressed by its midpoint and radius. For example, if m = mid(X), 
r = rad(X), then we can write the formula above as X = midrad(m, r). 

We can also define the arithmetic operations over intervals. Let X = \x, x], Y = 

[y,y]eI(R), 

- X + Y = [x + y, x + y] 

- X -Y =[x-y,x-y\ 

- X Y = [min(xy,xy,xy,Ty),max(xy^,xy,xy,xy)} 



1 See reference [15], Section 11, Automatic differentiation. 



- X/Y=\x,x]-[l/y,l/y], Q^Y 

A vector is called an interval vector if all its components are intervals. Interval 
matrix can be similarly defined. For interval vectors and interval matrices, the 
concepts such as midpoint, width, radius, etc, and the arithmetic operations are 
defined in components. 

Let / : R n — > R be a function, if there exists an interval map 

F : I(R n ) -> I(R) 

such that for all Xi € Xi(i = 1,2, ... , n), 

F([xi,xi], [X2,X 2 ],. ■ ■ , [x n ,x n ]) = f(xi,x 2 , ■ ■ .,x„) 

holds, then we call F an interval expand of /. 

We call F : I(R n ) — > I(R) an interval map with inclusive monotonicity if 
XCY implies F(X) C F(Y) for any given intervals X and Y. The definitions 
above can all be extended to the situations in I(R n ) — > I(R n ). And it is easy 
to prove that all the polynomial operations satisfy the inclusive monotonicity. 

Krawczyk operator Krawczyk operator plays a key role in the real root veri- 
fication of interval arithmetic. The main accomplishment comes from the work 
of Krawczyk and Moore. We only list some important results here. Complete 
proofs can be found in [20]. 

Suppose / : D C R n — > R n is continuous differentiable on D. Consider the 
equation 

/(*) = (1) 

Let /' be the Jacobi matrix of /, F and F be the interval expand of / and /' 
with inclusive monotonicity, respectively. For X € 1(F)) and any y 6 X, define 
Krawczyk operator as: 

K(y, X)=y- Yf(y) + (I - YF'(X))(X - y) (2) 

where Y is any n x n non-singular matrix. 

Especially, we assign y = mid(X), so Formula (2) becomes 

K{X) = mid(X) - F/(mid(X)) + (I - FF'(X))rad(X)[-l, 1]. (3) 

Formula (3) is often used in practise. 

The reason why Krawczyk operator is so important is that it has some nice 
properties. 

Proposition 1 Suppose K(y,X) is defined as Formula (2), then 

1. If x* £ X is a root of Equation (1), then for any y £ X, we have x* 6 
K(y,X); 

2. For any y € X , if X n K(y, X) = holds, then there is no roots in X ; 



3. For any y £ X and any non-singular matrix Y , if K{y, X) C X holds, then 

Equation (1) has a solution in X ; 
4- Moreover, for any y £ X and any non-singular matrix Y , if K(y, X) is 

strict inclusive in X, then Equation (1) has only one root in X. 

With the properties above, we can easily develop a real root verification 
method, which is a little different from the classical one, and will be explained 
later in this paper. 

Meanwhile, with the hypothesis we set in introduction, all the systems con- 
sidered here are non-singular ones with only simple roots. So the Jacobian matrix 
at the zeros are all invertible. Thus, we often set Y = (mid.F'(-X")) _1 and the 
Krawczyk operator becomes 



This is also called the Moore form of Krawczyk operator. 
3 Real root isolation algorithm 

In this section, we will present our new algorithm for real root isolation based 
on hybrid computation. As mentioned before, our idea is to construct the initial 
intervals corresponding to the approximate zeros obtained by homotopy contin- 
uation, then verify them via Krawczyk interval iteration to obtain the isolation 
results. In the end, we combine these sub-procedures to give the final algorithm 
description. 

3.1 Construction of initial intervals 

To apply Krawczyk interval iteration, obviously the construction of initial inter- 
vals is a key procedure. We should guarantee both the correctness and efficiency, 
that is, make sure the initial box contains the corresponding accurate real root, 
and keep the interval radius as small as possible so as to shorten the iteration 
time. 

Thus a valid error estimation for the initial approximate zeros should be 
established. And we discuss this issue in both theory and practice aspects here. 

Error estimation theory The core problem of the construction of initial box 
is the choice of interval radius, which is indeed an error estimation for the ap- 
proximate zero. There are dozen of error analysis for this question, from classic 
results to modern ones, especially about the Newton method. For example, in 
[2], S. Smale et al. gave a detailed analysis. However, their method requires com- 
putation of high order derivatives, which is not so convenient for our problem. 
Here we employ Kantorovich Theorem to give our error estimation. 



K(X) 



mid(X) - (midF'(X))- 1 /(mid(X)) 

+ {I - (midF'(X))- 1 F'(X))rad(X)[-l, 1 



(4) 



Theorem 2 (Kantorovich) Let X and Y be Banach spaces, F : D C X — > Y 
be an operator, which is Frechet differ entiable on an open convex set Dq C D. 
For equation F{x) = 0, if the given approximate zero xq € Dq meets the following 
three conditions: 

1. ^'(xo) -1 exists, and there are real numbers B and r\ such that 

im^r 1 !! <b, hf'^o)- 1 ^)!! < v, 

2. F satisfies the Lipschitz condition on Dq: 

\\F'(x)-F'(y)\\ <K\\x-y\\, \fx,y e D , 



3. h = BKr] < i O(x , ^g^ T?) c D 0) 



2 

then we claim that: 



1. F(x) = has a root x* in 0(xq, 1 ~^l~ 2h -ri) C Dq, and the sequence {xk '■ 
Xk+i = Xk — F' (xk)^ 1 F(xk)} of Newton method converges to x* ; 

2. For the convergence of x* , we have: 



9(l-9 2k+1 ) 



Xk+i\\ < M , \. 2k+1 ' v (5) 



whe 



3. The root x* is unique in Dq H 0(xq, 1+v/ ^ 2h -v)- 

In the theorem, 0(x, r) denotes the ball neighborhood whose center is x and 
radius is r, and 0(x, r) refers to the closure of the ball neighborhood. The proof 
can be found in [8]. 

Since the approximation xq is already the result of homotopy process, what 
we care about is the initial interval w.r.t. xq, i.e. the proper upper bound for 
||x* — xo||. So we have the following proposition, which is a direct corollary of 
Kantorovich Theorem. 

Proposition 3 Let F = (/i, f%, ■ ■ . , f n ) T be a polynomial system, where f t G 
H[xi,X2, ■ ■ ■ ,x n ]. Denote by J the Jacobian matrix of F. For an approximation 
Xq £ K" , if the following conditions hold: 

1. J~ 1 (xq) exists, and there are real numbers B and rj such that 

IIJ^MI <B, \\J- 1 (x )F(x Q )\\ <v, 

2. There exists a ball neighbourhood 0(xq,lu) such that J(x) satisfies the Lip- 
schitz condition on it: 



\\J(x)-J(y)\\<K\\x-y\\, Vx.yeO^,.) 



3. Let h = BKr], 



1 , 1 - y/1 - 2h 

h < — , and ui > 77, 

~ 2 ~ h 



then F(x) = has only one root x* in 0(xq 1 lo). 

Proof. We consider F as an operator on R n R n , obviously it is Frechet 
diffcrcntiablc, and from 

F(x + h) = F{x) + J(x)h + o(h) 

we can get 

Hm \\F(x + h)-F(x)-J(x)h\\ = Q 

h^>0 \\h\\ 

Thus the first order Frechet derivative of F is just the Jacobian matrix J, i.e. 
F (x) = J(x). So by Theorem 2, the proof is completed immediately after 
checking the situation of \\x* — x \\. 



It is easy to know that 1 ~ v '^~ 2h < 2. So we can just assign lo = 2rj. Then we 
need to check whether K < in the neighborhood O(a;o,2?;). Even though 
the initial Xq does not satisfy the conditions, we can still find a proper x^ after 
several Newton iterations, since rj will approach zero, B is fixed and we only 
need to find an upper bound for the Lipschiz constant K. 



Constructive algorithm Now we will give a constructive procedure for the 
Lipschitz constant K. 

Let Jij = dfi/dxj, apply mean value theorem to each clement of J on 
O(x ,co) to get 

Jij(y) = Jij(x) + VJij(Cij) ■ (y - x), Vx,y e O(x ,uj) (6) 

where V = (d/dx\,d/dx2, ■ ■ ■ ,d/dx n ) is the gradient operator, and Cy is some 
point between x and y, i.e. Cij = x + K ij(y — x ) f° r some Kij £ [0, 1]. Setting 
(vJij(fe) • (y - x )) nxn = AJ > we have 

||AJ|| O0 = ||(|VJy(Ci i ).(|/- a ;)|) nXB || 0o 

< ll(l|VJ„-(Cy)ll2l|y-a;||2)„ x „lloo 

< n||(||VJij(Ci3)llo ) TlXn l|cx> • \\y - x \\oc 

n 

= n max V ||VJij(C»j)||oo ■ \\y - ^IU- (7) 

±<%<n * — * 

" ~ i=i 

Note that VJij(Qj) is a vector, so if we use | • | max to denote the maximum 
module component of a vector, then we have 

n 

II AJ||oo < n ■ max N |VJy (C*j ) I max ' || V «E||oo* 

(8) 

l<t<n c — * 
" " 3=1 



Let Hi = ( dx-dxk )" x " ^ c ^ e Hessian matrix of /j, and let Hi = (h[ l \ . . . , h$), 
where hj are the column vectors. Then we have 

n 

HAJHoo^n max V \hf{Cij)U™ • \\v ~ x\\oo- (9) 

Ki<n L — * J 
~ ~ 3=1 

For convenience, we construct Xq = midrad(a;o, oS) with Xq as centre and 
lo = 1r\ as radius. 

Now we have \hf' (dj)\ max < \hf\x )\ max . So 

n 

IIAJHoo <n max V |/if(X )| max ' || 2/ 33 1 1 oo • 

(10) 

Ki<n — ' J 

Therefore 

n 

X = n max V \hf (X )\ niax (11) 

~ ~ 0=1 

is the Lipschitz constant w.r.t. J. 

Now we give an algorithm for computing initial intervals in Algorithm 1. 



Algorithm 1 init_width 

Input: Equation F; Approximation xq; Number of variables n 
Output: Initial interval's radius r 
1: repeat 

2: x = x - J~ 1 (x )F(x ); 

3: B = || J-^sboJIIwJ V = \\J- 1 (x )F(x )\\ oo ; 

4: uj = 2rj- 

5: X[) = midrad(cco, w); 

6: # = 0; 

7: for i = 1 to n do 

8: Compute the Hessian matrix Hi — (h\ , h\ , ■ ■ ■ , ft« ) of F on Xo; 

9: if E;=i \hf (Xo) | max > A" then 

10: K = J2%i\hf\X )\ m ^; 

11: end if 

12: end for 

13: K — n - K; 

14: ft = Btfjj; 

15: until ft < 1/2 

16: return r = 1 ~ v }~ 2h ri 



3.2 Empirical estimation 



As so far, we have established a rigorous method to construct initial intervals. 
This method takes a complex approximate zero as input to obtain an initial box. 
But in practice we often find many approximations with "large" imaginary parts 
which strongly indicate that they are non-real. A natural question is 

Can we detect these non-real roots without using interval arithmetic? 

Let z be an approximation of the real root Because 

||£Ke(z)-£|| < ||*- 

then we can see the real part 5Hc(z) is also an approximation of this root and is 
even closer. So we can simply replace Xq by D\c(xq) in Algorithm 1 to construct 
the initial box. 

The other consideration is the efficiency of numerical computation. When we 
use Proposition 3, lots of interval matrix operations would be executed, which 
cost much time than the point operations. So if we can find an empirical estimate 
radius, which can be computed much faster, but is still valid for most of the 
equations, then that will be a good choice in practice. 

We now give one such empirical estimation. 

For F — 0, let x* be an accurate root and Xq be its approximation. Although 
the mean value theorem is not valid in complex space, the Taylor expand is still 
valid. And the polynomial systems considered here are all continuous, so we 
suppose the equation satisfies the mean value theorem approximately: 

= F(x*)ttF(x o ) + J(S)(x*-x o ) (12) 

where £ is between x* and Xq. So we have 

x* -Xo^-J-'i^Fixo). 

Let J(£) = J(x ) + A J, then 

J(0 = J(x Q )(I + J- 1 (x Q )AJ), 
J- x (0 = (/ + J- 1 (x Q )AJ)- 1 J- 1 (x Q ). (13) 

For A J, we can get an estimation similar to Formula (9): 

71 

IIAJHoo < n max V \hf (Climax -\W - Moo- (14) 
Ki<n *• — * J 

From our hypothesis, x* and xq are very close, so are Qj and xq. Thus, we 
approximate Xq with Meanwhile, from Xq, after a Newton iteration, we get 
Xi = x — J~ 1 (x )F(x ). Thus we may consider that the distance between x* 
and Xo is more or less the same with that of Xq and x±, so we replace ||x* — Xq\\ 
with \\xi — x Q \\ = || J~ 1 (xq)F(x )\\ for approximation. 



So we get 

n 

IIAJHoo < n max V \hf(x )\ max ■ \\ J' 1 (x Q )F(x )\\ oc . (15) 

Ki<n — ' J 

~ ~ i=i 

Let A = maxi<i<„ YTj=\ Wj ( £C o)|max, then 

|| J- 1 (a; )AJ|| oo < nAHJ-^aoJIlLll-FC^Hoo. 

Because ||-F(a;o)||oo "C 1, the last formula is also far less than 1. So substitute 
that into Formula (13) we can get 

|| J-^Hoo < l|J " 1(:Eo)l100 . (16) 

Finally we obtain the empirical estimation 

||as* - ssolloo ~ ll^~ 1 (€)- F (*o)||oo 

- 1 -^A||J r - 1 ( ; ro)||?oll^(^o)||oo' 

Notice that the inequality (17) is only a non-rigorous estimation. All the com- 
putation in it are carried out in a point-wise way, so it is faster than Proposition 
3. In the numerical experiments later we will see that this empirical estimate 
radius performs very well. So we can use it to detect those non-real roots rather 
than the interval arithmetic. We describe that in Algorithm 2. 



Algorithm 2 iscomplex 

Input: Equation F; Approximation z; 

Output: true (z must be non-real), or false (z may be real). 
1: Compute Formula (17), denote the result by r'\ 
2: if any( 3m(z) > r' ) then 

3: return true; / / not a real root, continue to judge others 
4: else 

5: return false; // may be a real root, call interval arithmetic to verify 
6: end if 



In Algorithm 2, any() is a default function in Matlab, which returns true if 
there is non-zero component in a vector. 

3.3 Krawczyk-Moore interval iteration 

We now discuss about the real root verification with a given interval. In sec- 
tion 2.2, we have introduced the Krawczyk operator. With the properties in 



Proposition 1, we can determine whether an interval contains a real root by the 
relationship of the original interval and the one after the Krawczyk iteration. 

However, in practice, we can't expect the intervals to be entire inclusion or 
disjoint after just one iteration. Partly intersection is the most common cases 
that we encounter. Since the real root is still in the interval after the Krawczyk 
iteration, a normal method is to let X n K{X) be the new iteration interval. 
So suppose X (0) is the initial interval, the iteration rule is = n 

K{X( k) ), where K{X {k) ) is defined by Formula (4). This update rule can make 
sure that the size of X^ is non-increasing. But a problem is once we encounter 
K(X^) n X (k) == X {k \ the iteration will be trapped into endless loop. So we 
have to divide X^ k ' if this happened. 

Thus, we introduce a bisection function divide(). To ensure the convergence 
of our algorithm, we divide the longest dimension of an interval vector. 

This strategy may not be the optimal choice when the system's dimension is 
high. Greedy method or optimization algorithm will be studied in future work. 

We now give a formal description of divide function in Algorithm 3 and the 
Krawczyk-Moore iteration process in Algorithm 4. 



Algorithm 3 divide 
Input: Interval vector X 
Output: X w , X (2) 

1: Let Xi be the coordinate with the largest width in X 

2: X w = X; X (2) = X; 

3: X^ = Etof(X<) f mid(JCi)]; 

4: X (2) = [midfJCi),sup(X ( )]; 

5: return X (1) , X {2) 



3.4 Verification and refinement 

After the Krawczyk iteration, we already have all the real root isolated inter- 
vals, but these are not the final results. Since we require an isolation of disjoint 
intervals, we have to check the possible overlaps. 

On the other hand, some intervals may not be as small as required by users, 
so we can narrow them via bisection method until they match the requirement. 

We discuss these details in this subsection. 



Remove the overlaps There is a basic hypothesis: for non-singular systems, 
each root has an approximation, and from this approximation, the iteration will 
end up in its corresponding accurate root, not the other zeros. So we only have 
to remove the overlaps, and the number of real roots won't change. 

However, we want to expand our algorithm into multi-roots cases. And in 
that situation, it is possible that two isolated intervals contain the same real 



Algorithm 4 Krawczyk 



Input: F; initial box X; isolation boxes real jroots; number of real roots ureal 
Output: symbol of whether there is a real root flag; realjroots; ureal 

1: Y = midC-F'CX))- 1 ; X t = K(X), where K(X) is define by Formula (4); 

2: if X t n X == then 

3: return flag = false; 

4: end if 

5: while not (X t C X) do 
6: if X t n X == X then 

7: [X (1) ,X (2) ] = divide(X); 

8: [fl, realjroots, ureal] = Krawczyk(F, X^ , realjroots, ureal); 

9: if /l == false then 

10: [/2, real jroots, ureal] — Krawczyk(_F, X' 2 ' , real jroots, ureal); 

11: end if 

12: return fl or /2; 

13: end if 

14: X = X t C\X; 

15: Y" = (midF'(X))- 1 . 

16: X t = K{X); 

17: if X t n X == then 

18: return flag = false; 

19: end if 

20: end while 

21: ureal = ureal + 1; 

22: realjroots[nreal] = Xt 

23: return flag = true, reaZ jroots, ureal; 



zero. So whether or not the overlap part contains a real root, our algorithm has 
its corresponding processes. See Algorithm 5 for details. 

The function KrawczykQ in Algorithm 5 is a little bit different from that in 
Krawczyk-Moore iteration. In Krawczyk-Moore iteration, we have to store the 
information of isolated real root intervals, so the real jroots and ureal are in the 
function arguments. However, we only need to know whether there is a real root 
here, so only the symbol variable flag is returned. The situation is the same in 
Algorithm 6. 

Narrow the width of intervals As we said in the introduction, the real 
root isolation can be viewed as a kind of solving equations. And the width 
of the isolated intervals is just like the accuracy of solutions. Similar to the 
former algorithms, we can require the program to return an answer in specified 
range. The difference is the symbolic algorithm can get any precision they want 
in theory, but our floating point number calculation can't beat the machine 
precision. In fact, in the Matlab environment that we implement our algorithm, 
the results width won't be smaller than the system zero threshold 1 . 



In Matlab2008b that we do the experiments, the zero threshold is 2.2204e-016. 



Algorithm 5 disjoint_process 

Input: Isolated intervals realjroots; number of real roots ureal; F 
Output: Checked isolated intervals realjroots; ureal 
1: k = 0; 

2: for i = 1 to ureal do 

3: = reaLroots [i] ; newjroot = true; 

4: for j = 1 to k do 

5: y = real_roots[j]; 

6: z = xny ; 

7: if Z == then 

8: continue; 

9: end if 

10: flag = Krawczyk(F,Z); 

11: if flag == true then 

12: newjroot = false ; break; 

13: else 

14: X = X \ Z; 

15: realjroats\j] = reaLroois [j] \ .Z; 

16: end if 

17: end for 

18: if newjroot == true then 

19: k — k + 1 ; reaLroois [fc] = X 

20: end if 

21: end for 

22: return realjroots ,nreal = k; 



Wc also use bisection to do the narrowing job. Since there is only one root 
in the interval, we only have to continue dividing and checking the half that 
contains that zero. Formal description is in Algorithm 6. 

3.5 Algorithm description 

Up to now, we have discussed all the parts of real root isolation algorithm in 
detail. We give the final main program in Algorithm 7. 

4 Experiments 

Now we apply our new method to some polynomial systems and do some com- 
parison with some former algorithms. 

All the experiments are undertaken in Matlab2008b, with Intlab [14] of Ver- 
sion 6. We use Hom4ps-2.0 [12] as our homotopy continuation tool to obtain 
initial approximate zeros. Since computation time will be listed below, the com- 
puter information is also given here: OS: Windows Vista, CPU: Inter®Core 2 
Duo T6500 2.10GHz, Memory: 2G. 



Algorithm 6 narrowing 



Input: Isolated intervals realjroots; Number of real roots ureal; F; Threshold r 
Output: realjroots 



1 


for i = 1 to ureal do 


2 


X — reaLroots[i]; 


3 


while any(rad(X) > r) do 


4 


[Y (1) ,Y (2) ] = divide(X 


5 


flag = Krawczyk(F, Y (1) 


6 


if flag == true then 


7 


X = Y^; 


8 


else 


9 


X = Y^ 


10 


end if 


11 


end while 


12 


real_roots[i] — X; 


13 


end for 


14 


return reaLroots 



4.1 Demo example 

We begin our illustration by a simple example. 

Example 1 Consider the real root isolation of the system below. 

x 3 y 2 + x + 3 = 
Ayz 5 + %x 2 y i z A -1 = 
x + y + z — 1 = 

The homotopy program tells us this system has 28 complex zeros in total. And 
we get the following results after calling our real_root_isolate program, 
intval = 

[ -0.94561016957416, -0.94561016957415] 
[ 1.55873837303161, 1.55873837303162] 
[ 0.38687179654254, 0.38687179654255] 
intval = 

[ -1.18134319868123, -1.18134319868122] 
[ -1.05029487815439, -1.05029487815438] 
[ 3.23163807683560, 3.23163807683561] 
intval = 

[ -2.99999838968782, -2.99999838968781] 
[ 0.00024421565895, 0.00024421565896] 
[ 3.99975417402886, 3.99975417402887] 
intval = 

[ -0.79151164911096, -0.79151164911095] 
[ 2.11038450699949, 2.11038450699950] 
[ -0.31887285788855, -0.31887285788854] 
The order of variables : 



Algorithm 7 rcal_root Jsolate 

Input: Equation F(x); number of variables n; Threshold r; 

Output: Isolated intervals of F(x) = and number of real roots ureal 

1: [complex jroots,ncomplex] = hom4ps(J ? ,n); 

2: Initialize real-roots to be empty; ureal = 0; 

3: for i = 1 to ncomplex do 

4: z = complex jroots[i]; 

5: if iscomplex(F, z) then 

6: continue; 

7: end if 

8: r — init_width[F,z,n]; 
9: X = midrad(SHe(z), r); 

10: [flag, real-roots, ureal] = Krawczyk(F, Xo, realjroots, ureal); 
11: end for 

12: [realjroots, ureal] = disjoint_process(reaLroots, ureal, F ); 
13: realjroots = narrowing(reaLro(rfs, ureal, F, r); 
14: return realjroots, ureal; 



'y' 

'z' 

The number of real roots : 4 

We verify the answers above with the DISCOVERER [22] package under 
Maple, which also return 4 isolated real roots. Here we show its output in float- 
ing point number format, i.e. 

[[-2.999998391, -2.999998389], [0.2442132427e-3, 0.2442180230e-3], [3.999754090, 3.999754249]], 
[[-1.181343199, -1.181343199], [-1.050294975, -1.050294818], [3.231637836, 3.231638372]], 
[[-.9456101805, -.9456101656], [1.558738033, 1.558738728], [.3868716359, .3868719935]], 
[[-.7915116549, -.7915116400], [2.110384024, 2.110385000], [-.3188729882, -.3188727498]]. 
And we can see the answers perfectly match the ones of our program. 

We list some information during the calculation of our algorithm here for 
reference. Only the 4 real ones are given, and the other non-real ones are all 
detected by our empirical estimation method. We mention that all the imagi- 
nary parts of complex roots are significant larger than the initial radius of our 
algorithm in order of magnitude in this example. 





rootl 


root2 


root3 


root4 


B 


1.060227 


1.192159 


2.000864 


0.874354 


K 


14.941946 


7.198937c+003 


4.095991e+003 


16.988990 


V 


4.260422e-016 


4.20807e-016 


8.882333e-016 


5.764449e-016 


h 


2.024791e-014 


1.083446e-011 


2.183861e-011 


2.568823e-014 


estimate-rad 


4.274976e-016 


4.208067e-016 


8.882344e-016 


5.779921e-016 


cmpirical-rad 


1.015249e-015 


1.29164e-012 


2.156138e-015 


1.559270e-015 



Table 1. Key quantities comparison 



We give some remarks on table 1. In the first row, rootl to root4 are refer to 
the 4 real roots mentioned above respectively. And B, K, 77, h are exactly the same 
as they are defined in algorithm 1. The estimate-rad are the radius obtained via 
algorithm 1, while the empirical-rad are refer to the ones calculated by Formula 
(17). 

We say a little more words about the empirical-rad. Firstly, although the 
empirical ones are basically larger than the rigorous error radius, they are still 
small enough, which hardly have any influence on the efficiency of interval it- 
eration. We will see this in the comparison experiments later. But avoiding of 
interval matrix computation is very helpful to the algorithm. Secondly, the ra- 
dius obtained from algorithm 1 are so small that they are even comparable to 
the zero threshold of Matlab system 1 . And this could bring some uncertainty 
of floating point operation to our algorithm, such as misjudgement of interval 
inclusion in Intlab, etc. So we intend to use empirical estimation bound in next 
experiments. 

For cyclic6 (see Appendix A), the classic symbolic algorithm can do nothing 
due to the difficulty of triangularization. Meanwhile, we can easily get the 24 
isolated real roots intervals with our real_root_isolate program. 

4.2 Comparison experiment 

Many benchmarks have been checked with our realjroot_isolate program. 
Since the time complexity of both triangularization and homotopy continuation 
are difficult to be analyzed, we mainly focus on the isolation results and the 
program execution time. 

We investigate over 130 benchmarks provided by Hom4ps [1], among which 
about 40 equations are non-singular systems. We apply our program to these 
equations and all the experiments receive the right answers. Here we list a few 
of them (see Appendix A for details). 

The column real roots in Table 2 tells the number of intervals that our pro- 
gram isolated. Compared with the results of DISCOVERER, the new algorithm 
indeed works out all equations that are beyond the capability of classic symbolic 
algorithm. Moreover, the last column show that our empirical estimate method 
detects all the non-real roots successfully. 

Table 3 gives the comparison of program execution time. The total time is not 
equal to the sum of homotopy time and interval iteration time because we only 
count the CPU time, and there are other tasks such as I/O, format transform, 
etc. 

Table 3 also shows that interval iterations consume more time than homotopy 
continuation. The reason is complicated and we enumerate some here: 

1. The homotopy focus only on floating point number, while Krawczyk cares 
about intervals; 

1 As mentioned before, the zero threshold in Matfab2008b is 2.2204e-016, which is 
almost the same order of magnitude of those radiuses. 



Equation 


total roots 


real roots 


DISCOVERER 


complex roots detected 


barry 


20 


2 


2 


18 


cyclic5 


70 


10 


10 


60 


cyclic6 


156 


24 


N/A 


132 


desl8_3 


46 


6 


N/A 


40 


eco7 


32 


8 


8 


24 


eco8 


64 


8 


N/A 


56 


geneig 


10 


10 


N/A 





kinema 


40 


8 


N/A 


32 


reimer4 


36 


8 


8 


28 


reimer5 


144 


24 


N/A 


120 


virasoro 


256 


224 


N/A 


32 



Table 2. Real root isolation results comparison 



Equations 


Total time 


Homotopy time 


Interval time 


DISCOVERER 


barry 


0.421203 


0.093601 


0.327602 


0.063 


cyclic5 


2.948419 


0.218401 


2.652017 


0.624 


cyclic6 


9.984064 


0.639604 


9.063658 


N/A 


desl8_3 


4.180827 


0.702004 


3.385222 


N/A 


eco7 


2.371215 


0.265202 


2.012413 


15.881 


eco8 


3.946825 


0.499203 


3.354022 


N/A 


geneig 


4.243227 


0.249602 


3.868825 


N/A 


kinema 


3.946825 


1.014006 


2.808018 


N/A 


reimer4 


2.480416 


0.374402 


2.059213 


24.711 


reimer5 


12.963683 


3.073220 


9.578461 


N/A 


virasoro 


137.124879 


4.570829 


109.996305 


N/A 



Table 3. Execution time comparison, unit:s 



2. Hom4ps-2.0 is a software complied from language C, which is much more 
efficient than the tool that we use to implement our algorithm, say Matlab. 

3. The interval iteration time increases as roots number grows since we examine 
the approximate zeros one by one. So the parallel computation of homotopy 
is much faster. 

We believe that with efficient language such as C/C++, and parallel computa- 
tion, the implementation of our algorithm will be much faster. 

In order to verify our idea and see whether parallelization could help, we go 
into every approximate zero's iteration process. Some critical data are recorded 
in table 4. The avg. rad. of ans is the average radius of the final isolated intervals, 
while the avg. rad. of init. indicates the average radius of the initial intervals. 
The average time of each zero's interval iteration is shown in column avg. time 
of iteration along with the max interval iteration time in max time of iter. We 
think the consumption for each zero's process is acceptable. 

From table 4 we can see that the initial interval radiuses are extremely small, 
which leads to a nice process time for each iteration. We point out that almost 



Equation 


avg.rad.of ans 


avg.rad.of init. 


avg.time of iter 


max time of iter 


barry 


3.552714e-015 


1.377800e-014 


0.054600 


0.062400 


cyclic5 


1.614703e-009 


7.142857c-007 


0.113881 


0.140401 


cyclic6 


4.440892e-016 


2.137195e-015 


0.183951 


0.234002 


desl8_3 


3.768247e-007 


9.737288e-007 


0.241802 


0.296402 


eco7 


1.998401e-015 


1.483754e-013 


0.122851 


0.156001 


eco8 


2.109424e-015 


3.283379e-013 


0.183301 


0.218401 


geneig 


2.664535e-016 


5.721530e-014 


0.315122 


0.436803 


kinema 


1.998401e-015 


6.784427e-011 


0.157951 


0.218401 


reimer4 


1.110223e-016 


1.258465e-014 


0.122851 


0.156001 


reimer5 


1.110223e-016 


4.754080e-014 


0.195001 


0.421203 


virasoro 


9.472120e-009 


2.265625e-006 


0.387844 


0.624004 



Table 4. Detail data for each iteration, unit:s 



all real root checks are done by just one Krawczyk iteration, and hardly any 
overlap is found after all the Krawczyk iteration processes due to the small 
initial intervals that we give. All of these save significant execute time of our 
program. 

5 Conclusion 

For the non-singular polynomial systems with variables' number equals equa- 
tions' number, this paper presents a new algorithm for real root isolation based 
on hybrid computation. The algorithm first applies homotopy continuation to 
obtain all the initial approximate zeros of the system. For each approximate zero, 
an initial interval which contains the corresponding accurate root is constructed. 
Then Krawczyk operator is called to verify all the initial intervals so as to get all 
the real root isolated boxes. Some necessary check and refinement work are done 
after that to ensure the boxes are pairwise disjoint and meet width requirement. 

In the construction of initial intervals, we give a rigorous radius error bound 
based on the corollary of Kantorovich theorem. Some constructive algorithms 
are presented for both real and complex approximate zeros. Meanwhile, we in- 
troduce an empirical estimate radius, which has a nice performance in numerical 
experiments. 

In the modification and implementation of Krawczyk iteration algorithm, 
some problems of interval arithmetic are also discussed in this paper. 

At last we utilize some existing tools to implement our algorithm under Mat- 
lab enviornment. Many benchmarks have been checked along with compairson 
and analysis. 

We also mention some possible future work here. 

The construction of initial intervals is still too complicated and further opti- 
mization shall be studied. Also the empirical estimation with more efficiency and 
accuracy is a considerable question. The divide strategy in Krawczyk iteration 
could also be improved, which may be helpful in the high dimension cases. 



In the aspect of implementation, replacing the Matlab implementation with 
C/CH — h codes may improve the performance of our program. Parallel computa- 
tion can solve another bottleneck of our problem. And for some small systems, 
or equations with special property, the classic symbolic algorithm could be even 
faster. So the tradeoff of symbolic and numerical computation is also an inter- 
esting direction. 
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A Benchmarks 

1. barry: Number of variables:3, Number of equations:3,Max degree:5 

-x 5 + y 5 - 3y - 1 = 
5y 4 - 3 = 
-20x + y- z = Q 

2. cyclic5: Number of variables:5, Number of equations:5,Max degree:5 

Xi + X 2 + X 3 + X4 + X 5 = 
X1X2 + X 2 X 3 + X3X4 + X4X5 + XiX 5 = 
X1X2X3 + X2X3X4 + X3X4X5 + X4X5X1 + X5X1X2 = 
X1X2X3X4 + X2X3X4X5 + 2:3X4X52:1 + X4X5X1X2 + X5X1X2X3 = 

X1X2X3X4X5 — 1 = 



3. cyclic6: Number of variables:6, Number of equations:6,Max degree:6 

.T + Xi + X 2 + X 3 + X 4 + X 5 = 
XqXi + X\X 2 + X2X3 + X3X4 + X4X5 + X$Xq = 
XqX\X2 + X1X2X3 + X2X3X4 + £3X4X5 + X4X5X0 + X5X0X1 = 
X XiX 2 X 3 + X1X2X3X4 + X2X3X4X5 + X3X4X5X0 

+X4X5X0X1 + X5X0X1X2 = 
X0X1X2X3X4 + X1X2X3X4X5 + X2X3X4X5X0 + X3X4X5X0X1 

+X4X5X0X1X2 + X5X0X1X2X3 = 
X0X1X2X3X4X5 — 1 = 

4. desl8_3: Number of variables:8, Number of equations:8,Max degree:3 

15a 33 ai a2i - 162a 10 a 2 2 - 312ai a 2 o + 24ai a 30 + 27a3ia 2 i 

+24a 32 a 2 o + 18022010032 + 30a 2 2O30 + 84a 3 iai = 
28a 22 aioa 3 3 + 192a 30 + 128a 32 a 10 + 36a 31 a 20 + 36a 33 a 2 o 

— 300ai a 2 i + 40a 32 a2i — 648a^ + 44a 2 2a 3 i = 
180a 33 ai - 284a 2 2aio - 162ai + 60022032 + 50a 32 aio 

+70a 30 + 55a 33 a2i + 260a 3 i - 112a 2 o = 
6a 3 3ai a 2 o + 10a 22 ai a3i + 8a 32 aioa2i - 162oi a 2 i 

+16a 2 ia 30 + 14a 3 ia 2 o + 48a w a 30 = 
4a 2 2aioa 3 o + 2a 3 2ai a 2 o + 6a 2 oa 30 - 162aj a 2 o + 3a 3 ia 2 iai = 
66033010 + 336032 + 90a3i + 78022033 — 1056aio — 90a2i = 
-240aio + 420a 33 - 64a 22 + H2a 32 = 
136a 33 - 136 = 

5. eco7: Number of variables:7, Number of equations:7,Max degree:3 

X7X1 + X7X1X2 + X7X2X3 + X7X3X4 + X7X4X5 + X7X5X6 — 1=0 
X7X2 + X7X1X3 + X7X2X4 + X7X3X5 + X7X6X4 — 2 = 
X7X3 + X7X1X4 + X7X2X5 + X7X6X3 — 3 = 
X7X4 + X7X1X5 + X7X2X6 — 4 = 
X7X5 + X7X1X6 — 5 = 
X6X7 — 6 = 
xi + X2 + X3 + X4 + X5 + + 1 = 



6. eco8: Number of variables:8, Number of equations:8,Max degree:3 

XiX$ + X 8 XiX 2 + X 8 X 2 X 3 + X 8 X 3 X4 + X8X4X5 

+x 8 x 5 x 6 + x 8 x 6 x 7 -1 = 

X 2 X 8 + X S XiX 3 + X S X 2 X4 + X8X3X5 + X 8 X 6 X4 + X S X7X 5 -2 = 
X8X3 + X S XiX4 + X 8 X 2 X 5 + X 8 X 6 X 3 + X8X7X4 — 3 = 
X S X4 + X 8 XiX 5 + X 8 X 2 X 6 + X S X 7 X 3 —4 = 

x 8 x 5 + x 8 xix e + x 8 x 7 x 2 -5 = 
x^xq + x$x 7 xi —6 = 
x 7 xg, — 7 = 

%i + x 2 + x 3 + X4 + x 5 + x 6 + x 7 + 1 = 

7. geneig: Number of variables:6, Number of equations: 6, Max degree:3 

— lOziXg + 2x 2 Xg — x 3 x\ + X4x\ + 3x 5 Xg + xix 6 + 2x 2 x 6 

+x 3 x 6 + 2x4X 6 + x 5 x 6 + 10xi + 2x 2 — x 3 + 2x 4 — 2x 5 = 
2xixl — \\x 2 x\ + 2x 3 x\ — 2x4x1 + x 5 x% + 2xix 6 + x 2 x 6 

+2x 3 Xq + X4X 6 + 3x 5 x 6 + 2x\ + 9x 2 + 3x 3 — X4 — 2x$ = 
— XiXq + 2x 2 Xq — \2x 3 x\ — X4x\ + x$Xq + xiXq + 2x 2 Xq 

— 2x4Xq — 2X5X6 — X\ + 3X2 + IOX3 + 2x4 — £5 = 

xiXg — 2x 2 Xg — x 3 Xg — 10x4Xg + 2x 5 Xg + 2xix 6 + x 2 x 6 
— 2X3X6 + 2X4X6 + 3X5X6 + 2xi — x 2 + 2x3 + 12x4 + X5 = 
3xiXg + x 2 Xg + x 3 Xg + 2x 4 Xg — llx5Xg + X1X6 + 3x 2 Xe 
—2x3X6 + 3x4X6 + 3x5X6 — 2xi — 2x2 ^ x 3 + X4 + 10x5 = 

xi + x 2 + X3 + X4 + X5 — 1 = 

8. kinema: Number of variables:9, Number of equations:!), Max degree:2 

z\ + zl + z 3 - 12zi - 68 = 
z\ + z\ + z\ - 12z 5 - 68 = 
Zj + zl + zl - 24z 8 - 12z 9 + 100 = 
Z1Z4 + z 2 z 5 + z 3 z 6 — 6zi — 6z 5 — 52 = 
Z1Z7 + z 2 z% + z 3 zq — 6zi — 12zg — 6zq + 64 = 
Z4Z7 + z 5 z s + z 6 z 9 — 625 — 12z 8 — 6zg + 32 = 
2^2 + 2^3 — Z4 — z$ — 2zq — z 7 — zq + 18 = 
zi + z 2 + 2z 3 + 2z4 + 2zq — 2z 7 + zg — zg — 38 = 
z\ + z 3 - 2z4 + z 5 — z 6 + 2z 7 — 2z$ + 8 = 



9. reimer4: Number of variables:4, Number of equations:4,Max degree: 5 

2x1 - 2x i + 2xl-2xj-l=0 
2x-^ ■ H - 2x^ 2x^ ■ 1 = 
2x^ 2^2 ~\~ 2x^ 2x^ ■ 1 = 
2x^ 2^2 "I - 2^3 — 2X4 ■ 1 — — 



10. reimer5: Number of variables:5, Number of equations:5,Max degree: 6 



2x\ 


2x 2 


Y2x\ 


- 2x 2 4 ~ 


-2x 2 5 - 


1 


= 


2x\ 


— 2x\ - 


V2x% 


-2x1- 


-2x1- 


1 


= 


2x\ 


— 2Xn ~ 


Y2x\ 


-2x\- 


-2x\- 


1 


= 


2x\ 


— 2x 2 - 


Y2x\ 


-2x1- 


-2x1- 


1 


= 


2x1 


— 2x\ - 


v2x% 


-2x1- 


-2x1- 


1 


= 



11. virasoro: Number of variables:8, Number of equations:8,Max degree:2 

2x\X4 — 2x\x-j + 2x 2 X4 — 2x2Xq + 2x3X4 — 2x3X5 + 8x 4 

+8x4X5 + 2x4X6 + 2x4X7 + 6x4X8 — 6x5X8 — X4 = 
2x4X5 — 2x\Xq + 2x2X5 — 2x2X7 — 2x3X4 + 2x3X5 + 8x4X5 

—6x4X8 + 8x5 + 2x5X6 + 2x5X7 + 6x5X8 — X5 = 
— 2X4X5 + 2xiX6 — 2X2X4 + 2X2X6 + 2x3X6 — 2x3X7 + 2X4X6 

+ 2X5X6 + 8X5 + 8X6X7 + 6X6X8 — 6x7X8 — X6 = 

— 2X1X4 + 2X1X7 — 2X2X5 + 2X2X7 — 2X3X6 + 2X3X7 + 2X4X7 

+ 2X5X7 + 8X6X7 — 6X6X8 + 8X7 + 6x7X8 — X7 = 
8x^ + 8x1X2 + 8x1X3 + 2xiX4 + 2xiX5 + 2xiX6 

+2x1X7 — 8x2X3 — 2x4x7 — 2x5X6 — xi = 
8X1X2 — 8X1X3 + 8X2 + 8X2X3 + 2X2X4 + 2X2X5 

+ 2X2X6 + 2X2X7 — 2X4X6 — 2X5X7 — X2 = 

—8X1X2 + 8X1X3 + 8X2X3 + 8X3 + 2X3X4 + 2X3X5 

+ 2X3X6 + 2X3X7 — 2X4X5 — 2X6X7 — X3 = 
— 6X4X5 + 6x4X8 + 6X5X8 — 6X6X7 + 6X6X8 + 6x7X8 + 8x| — X8 = 



